#Plot of core and edge proportions (plus sds) with no parasite cost at k=32000 and cell.size =5

setwd("/Users/ben/Documents/Papers/Submitted/Rhabdias/Parasite lag model/Pathol Freq Dep Dyn Load/Var R")


y.error.bars<-function(x, y, ybar, lty="solid", ...){
	points(x, y, ...)
	arrows(x, y-ybar, x, y+ybar, code=3, angle=90, length=0.1, lty=lty, ...)
}

data1<-read.table(file="Pathol Freq Depdt SI R1 Final.txt", header=T)
names(data1)[3]<-"beta"
attach(data1)

par(mfrow=c(2, 2)) #resize to 800 x 800

## First plot
core<-tapply(core.prop, list(surv.cost, beta), mean, na.rm=T)
edge<-tapply(edge.prop, list(surv.cost, beta), mean)
core.sd<-tapply(core.prop, list(surv.cost, beta), sd)
edge.sd<-tapply(edge.prop, list(surv.cost, beta), sd)
x<-as.numeric(colnames(core))

plot(c(min(x), max(x)), 0:1, ,type="n", cex.lab=1.5, cex.main=1.5, bty="l", xlab="", ylab="Observed prevalence", main="r=1")

cols<-c("black", "grey40", "grey70")
for (i in 1:3){
	if (i==1) y.error.bars(x, core[i,], core.sd[i,], lty="solid", pch=19,cex=1.5, col=cols[i])
	else points(x, core[i,], pch=19, col=cols[i], cex=1.5)
	lines(x, core[i,], col=cols[i])
}
for (i in 1:3){
	if (i==1) y.error.bars(x, edge[i,], edge.sd[i,], lty="dashed", pch=21,cex=1.5, col=cols[i])
	else points(x, edge[i,], pch=21, col=cols[i], cex=1.5)
	lines(x, edge[i,], col=cols[i])
}


legend(0.15, 0.6, legend=c("Range core", "Range edge"), pch=c(19, 21),bty="n")
legend(0.15, 0.8, legend=c("=0", "=0.1", "=0.4"), lty="solid",col=cols, bty="n", title="Cost of parasitism")

detach(data1)

data2<-read.table(file="Pathol Freq Depdt SI R3 Final.txt", header=T)
names(data2)[3]<-"beta"
data2<-data2[-1,]
attach(data2)

core<-tapply(core.prop, list(surv.cost, beta), mean, na.rm=T)
edge<-tapply(edge.prop, list(surv.cost, beta), mean)
core.sd<-tapply(core.prop, list(surv.cost, beta), sd)
edge.sd<-tapply(edge.prop, list(surv.cost, beta), sd)
x<-as.numeric(colnames(core))

plot(c(min(x), max(x)), 0:1, ,type="n", cex.lab=1.5, cex.main=1.5,  bty="l", xlab="", ylab="", main="r=3")

cols<-c("black", "grey40", "grey70")
for (i in 1:3){
	if (i==1) y.error.bars(x, core[i,], core.sd[i,], lty="solid", pch=19,cex=1.5, col=cols[i])
	else points(x, core[i,], pch=19, col=cols[i], cex=1.5)
	lines(x, core[i,], col=cols[i])
}
for (i in 1:3){
	if (i==1) y.error.bars(x, edge[i,], edge.sd[i,], lty="dashed", pch=21,cex=1.5, col=cols[i])
	else points(x, edge[i,], pch=21, col=cols[i], cex=1.5)
	lines(x, edge[i,], col=cols[i])
}


legend(0.15, 0.6, legend=c("Range core", "Range edge"), pch=c(19, 21),bty="n")
legend(0.15, 0.8, legend=c("=0", "=0.1", "=0.4"), lty="solid",col=cols, bty="n", title="Cost of parasitism")
detach(data2)

## Second plot

attach(data1)
summary<-tapply(x.diff, list(surv.cost, beta), mean)
summary.sd<-tapply(x.diff, list(surv.cost, beta), sd)
summary[which(summary==Inf)]<-NA
summary.sd[which(is.nan(summary.sd)==T)]<-NA
x<-as.numeric(colnames(summary))

plot(c(min(x), max(x)), c(-4, max(summary, na.rm=T)), ,type="n", cex.lab=1.5, cex.main=1.5,  bty="l", xlab="Transmission rate", ylab="Lag distance (arbitrary units)")
for(i in 1:3){
	if (i==1) y.error.bars(x, summary[i,], summary.sd[i,], lty="solid", col=cols[i], pch=21, cex=1.5)
	else points(x, summary[i,], col=cols[i], pch=21, cex=1.5)
	lines(x, summary[i,], col=cols[i])
}
legend("topright", legend=c("=0", "=0.1", "=0.4"), lty="solid",col=cols, bty="n", title="Cost of parasitism")
detach(data1)
attach(data2)
summary<-tapply(x.diff, list(surv.cost, beta), mean)
summary.sd<-tapply(x.diff, list(surv.cost, beta), sd)
summary[which(summary==Inf)]<-NA
summary.sd[which(is.nan(summary.sd)==T)]<-NA
x<-as.numeric(colnames(summary))

plot(c(min(x), max(x)), c(-4, max(summary, na.rm=T)), ,type="n",cex.lab=1.5, cex.main=1.5, bty="l", xlab="Transmission rate", ylab="")
for(i in 1:3){
	if (i==1) y.error.bars(x, summary[i,], summary.sd[i,], lty="solid", col=cols[i], pch=21, cex=1.5)
	else points(x, summary[i,], col=cols[i], pch=21, cex=1.5)
	lines(x, summary[i,], col=cols[i])
}
legend("topright", legend=c("=0", "=0.1", "=0.4"), lty="solid",col=cols, bty="n", title="Cost of parasitism")
detach(data2)

## Clean up

rm(list=ls())

